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Abstract 

Discrete particle simulation, a combined approach of computational fluid dynamics and discrete 
methods such as DEM (Discrete Element Method), SPH (Smoothed Particle Hydrodynamics), 
PIC (Particle-In-Cell), etc., is becoming a practical tool for exploring lab-scale gas-solid systems 
owing to the fast development of its parallel computation. However, the gas-solid coupling and 
the corresponding fluid flow solver remain immature. In this work, we presented a modified 
lattice Boltzmann approach to consider the effect of both the local solid volume fraction and the 
local relative velocity between the particles and the fluid, which was different from the traditional 
volume-averaged Navier-Stokes equations. This approach is combined with a time-driven hard 
sphere algorithm to simulate the motion of individual particles in which particles interact with 
each other via hard-sphere collisions but the collision detection and motion of the particle are 
performed at constant time intervals, and the EMMS (energy minimization multi-scale) drag has 
been coupled with lattice Boltzmann based discrete particle simulation to improve its accuracy. 
Two typical fluidization processes, namely injection of a single bubble at incipient fluidization 
and fast fluidization in a riser, are simulated with this approach, showing a good agreement 
with published correlations and experimental data. The capability of the model to capture more 
detailed and intrinsic characteristics of particle-fluid systems are demonstrated. The method can 
be used straightforwardly with other solid phase solvers also. 

Keywords: Discrete particle simulation, Lattice Boltzmann method, Computational fluid 
dynamics, Fluidization, Multiphase flow, Simulation 



1. Introduction 

Gas-solid fluidization systems are widely encountered in many industries for both physical 
and chemical processing. The earlier studies of these systems were focused on experimental 
investigations of their macroscopic hydrodynamic behavior to develop some corresponding cor- 
relations. With increasing computer capability, the computational fluid dynamics approach is 
adopted instead of experimental investigations in many cases, which has gained considerable 
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attention. In recent decades, significant advances have been made in hydrodynamic model- 
ing and simulation of gas-solid fluidization. Methods with varying scales and resolutions has 
been proposed, such as two-fluid model (TFM) (Anderson and Jackson, 1967; Ishii, 1975), dis- 
crete particle simulation (DPS) (Tsuji et al., 1993; Hoomans et al., 1996; Xu and Yu, 1997), 
quadrature-based moment methods (QBMM) (Desjardins et al., 2008; Fox, 2008), and direct 
numerical simulation (DNS) (Hu et al., 1992; Ma et al., 2006; Wang et al., 2010; Xiong et al., 

2012) . 

TFM treats the solid and the fluid as interpenetrating mixtures, each of which is governed 
by the conservation laws and requires constitutive equations for closure. In general, these con- 
stitutive equations are based on the kinetic theory of granular flow (Gidaspow, 1994) or simply 
empirical models, where large uncertainties still remain. Therefore, despite of the computational 
convenience and the wide application of TFM in multiphase flow simulation, its accuracy and ef- 
fectiveness are still in question in many circumstances. On the other hand, direct numerical simu- 
lation has been widely regarded as the most accurate method for the simulation of gas-solid flow. 
In this method, not only the motion of each individual solid particle is fully resolved, but also 
the hydrodynamic force acting on each individual solid particle is calculated directly from the 
stress on the fluid-solid boundary, together with the detailed flow field description around each 
particle. However, both spatial and temporal scales of DNS are limited by both the Kolmogorov 
microscale and the turbulence time scale for evolution, which leads to huge computational cost 
even at low Reynolds numbers. Therefore, though DNS is a useful tool for fundamental research, 
it is still formidable, at least in the foreseeable future, for predicting the hydrodynamics in large 
industrial scale fluidized beds. 

TFM is computationally more economic but inaccurate, while DNS is computationally more 
accurate but expensive. DPS, is somehow in between these two ends, and seems to give a good 
balance among accuracy, cost and efficiency. DPS has been proven to be effective in modeling 
various particle flow systems (Deen et al., 2007; Zhu et al., 2007, 2008), such as slugging flu- 
idized bed (Xu et al., 2007), spouted bed (Zhao et al., 2008), pneumatic conveying (Kuang et 
al., 2008), bubbling fluidized bed (Geng and Che, 2011), sound-assisted fluidized bed (Wang 
et al., 2011) and cyclone separator (Chu et al., 2011). In DPS, the continuum fluid flow is re- 
solved at the scale of computational cells in computational fluid dynamics (CFD), the motion 
of individual particles is described by the well-established Newton's equations of motion, and 
particle-particle interactions are presented by different collision models such as the hard-sphere 
model, the soft-sphere model. In DPS, fluid flows are traditionally described by the volume- 
averaged Navier-Stokes equations or their simplified forms (Tsuji et al., 1993; Hoomans et al., 
1996; Xu and Yu, 1997; Mikami et al., 1998), and those equations are commonly solved based 
on implicit schemes no matter by Fluent (Chu and Yu, 2008a, b; Chu et al., 2009a, b, 201 1; Wu et 
al., 2010; Zhao et al., 2010), OpenFOAM (Su et al., 201 1; Goniva et al., 2012), MFIX (Darabi et 
al., 201 1; Garg et al., 2012; Li and Guenther, 2012; Li et al., 2012a,b; Gopalakrishnan and Tafti, 

2013) , or in-house codes (Ouyang and Li, 1999a, b; Zhou et al., 2004a, b; Zhao et al., 2009; 
Wang et al., 2009; Wu et al., 2009). However, the parallel efficiency of those implicit solvers for 
fluid phase is not high, and they are hard to be implemented for the fast simulation of large-scale 
industrial systems. 

Lattice Boltzmann method (LBM) is an increasingly popular approach for solving complex 
flows due to its natural parallelism and algorithmic simplicity, etc. As an alternative compu- 
tational approach to the conventional CFD, LBM could be easily incorporated into DPS. Fil- 
ippova and Hanel (1997) combined lattice-BGK model and Lagrangian approach to simulate 
three-dimensional gas-particle flow through filters. However, it was based on one-way coupling, 
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that is, the effect of the fluid resistance on particles was considered without considering the cor- 
responding effect of particles on the fluid. Chen et al. (2004) simulated particle-laden flow over 
a backward-facing step with two-way coupling, where a modified lattice-BGK model was de- 
veloped for the fluid flow and a Lagrangian approach for particles, but they did not consider 
the effect of solid volume fraction on gas flows. Sungkorn et al. (2011) simulated turbulent 
gas-liquid bubbly flows with a relatively low gas holdup using stochastic Lagrangian model and 
LBM, which was a variant of DPS for gas-liquid flow to some extent. Specifically, they solved 
the continuous liquid phase by single-phase lattice Boltzmann equation (LBE) incorporated with 
large eddy simulation (LES) (Smagorinsky, 1963), excluding the gas volume fraction in the con- 
servation equations and its effect on drag force, and evolved the dispersed gas phase (i.e. the 
individual bubbles) by Lagrangian trajectories. 

In this paper, the equations of motion governing fluid flow in DPS are described by a modified 
LBE with a reasonable consideration of the effect of both the local solid volume fraction and the 
local relative velocity between particles and fluid, rather than the volume-averaged Navier-Stokes 
equations. The equations of motion governing individual particles are solved with time-driven 
hard-sphere (TDHS) model. It is noteworthy that the similar computational strategy has been 
implemented in DNS of particle-fluid systems (Wang et al., 2010) where the modified LBE was 
used with particles size much larger than the cell spacing, but in this paper, the partial saturation 
concept has been proposed to model objects much smaller than the cell spacing (i.e. porous 
media), and both the linear and the nonlinear drag effects of the solid phase (media) have been 
first considered so far in the lattice Boltzmann based discrete particle simulation. The proposed 
model is validated by two typical fluidization processes, namely a single bubble injected into a 
fluidized bed at incipient fluidization and particles transport in a circulating fluidized bed riser, 
with comparison to experimental data and published correlations, and demonstrating a good 
balance among accuracy, cost and efficiency. 



2. Numerical approach 

LBM, originally proposed by McNamara and Zantti (1988) as a smoothed alternative to lat- 
tice gas automata (LGA) , is an efficient second-order flow solver capable of solving various 
systems for hydrodynamics owing to its algorithmic simplicity, explicit solution of particle dis- 
tribution functions, natural parallelism and flexibility in boundary treatment. The objective of this 
research is to develop a lattice Boltzmann based discrete simulation approach for gas-solid sys- 
tems. For illustration, we used two-dimensional nine-velocity (D2Q9) lattice Boltzmann model 
as an example. As shown in Fig. 1, the governing equations of gas phase are solved based on 
D2Q9 lattice Boltzmann model and the solid particles distributed in the lattice cell are described 
by time-driven hard sphere model. 



2.1. LBM for gas phase hydrodynamics 

The particle distribution function / (f , x, v) is the number of particles per unit volume with 
velocity v at time t and position x, which is often used in non-equilibrium statistical mechanics. 
The evolution of the distribution function f is described by the Boltzmann equation (Chapman 
and Cowling, 1970): 

a/g,x,y )+Y a/(;,x,y )=Q 
dt dx 
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Figure 1: A schematic diagram of the D2Q9 computational fluid lattices and solid particles. 



The left-hand side of this equation is an advection term, while the right-hand side contains the 
collision operator O which describes the interaction between particles. External forces are ne- 
glected in this derivation. After a finite difference discretization of the Boltzmann equation in 
time, space and velocity, the LBE is simplified as (McNamara and Zanetti, 1988), 

fi(x + c,-M t + At) - fi(x, t) = (2) 

where the subscript i represents the orthogonal and diagonal directions in the Cartesian coordi- 
nates, fi is the particle distribution function in direction i, c, is the discrete velocity in direction i, 
x and t are the lattice simulation position and time, respectively, and £2, is collision term in direc- 
tion i. Bhatnagar, Gross and Krook developed a single relaxation time approximation, namely the 
so-called BGK model (Bhatnagar et al., 1954), in which the collision term £2, was approximated 
by 

Here r — 3v + y is the relaxation time related to kinematic viscosity and the discrete time step 
Af, jf q is the equilibrium distribution function defined as (Qian et al., 1992): 

\ / c,--u l(c--u) 2 lu-u\ , 

f, (p,u) = ^l + --- + --=- -----J, ,=0,1,2 ..■ 8 (4) 

where c, is defined as Co=0, and c,= {cos[(/- 1)^/2], sin[(/- 1)^/2]} for i= 1 ~ 4 and 
Cj= V2{cos [(/ - 5)n/2] , sin [(i — 5)n/2]} for i— 5 ~ 8. The weights are given by a»o= 4/9, 
a>;=l/9 for i= 1 ~ 4, djj-1/36 for i- 5 ~ 8, and c s is the speed of sound and equals to V3/3 in 
lattice unit, p and u are density and velocity, respectively. 

Through the Chapman-Enskog procedure, the incompressible Navier-Stokes equations can 
be obtained from the LBE in the limit of small Mach number (Sterling and Chen, 1996; Guo et 
al., 2002). The macroscopic quantities such as mass density and momentum density can then 
be obtained by evaluating the hydrodynamic moments of the distribution function fj(x, t) , and 
namely 

8 8 
P = £jfi and P u = Yj Ci f< (5) 

1=0 !=0 

The particle-fluid coupling is implemented by immersed moving boundary method (Noble 
and Torczynski, 1998; Cook et al., 2004), which introduces a term that depends on the percentage 
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of the cell saturated with fluid to modify the collision operator and represent the effect of both 
the solid volume fraction and the gas-solid slip velocity on hydrodynamics. The key idea of the 
immersed moving boundary method is that the effect of boundary replaced by the source term 
of force. Although most applications of the immersed moving boundary are used with objects 
much larger than the cell spacing (Noble and Torczynski, 1998; Cook et al., 2004; Feng et al., 
2007; Wang et al., 2010; Zhou et al., 2011; Xiong et al., 2012; Wei et al., 2013), the method is 
theoretically not confined to this scale. Here the partial saturation concept is extended to model 
objects much smaller than the cell spacing, and the LBE modified for partially saturated cells 
(Noble and Torczynski, 1998) reads, 

Mx + c,Af, t + At) = Mx, t) - -(1 - y(£s, t))C/Kx, - f^(x, t)) + y(e s , r)Q> (6) 

T 

where is an additional collision term that accounts for the bounce back of the non-equilibrium 
part of the distribution function. is given by (Noble and Torczynski, 1998) 

^ = /_,(x, t) - fi(x, t) + f*(p, v s ) - /!>, u) (7) 

where —i denotes the component of the distribution function opposite to i and v s is the average 
velocity of solid phase at the computational lattice. v s can be approximately calculated by particle 
number averaging as 

v. = 5£> (8) 

n lot 

where \ k is the velocity of solid particle k and « tot is the number of solid particles in fluid cell 
(see Fig.l). 

The weighting function y in Eq. (6) depends on the relaxation time r and the solid volume 
fraction s s of each cell, which is given by (Noble and Torczynski, 1998) 

Ss (t-0.5) ... 

T( £ s,t) = — — — (9) 

(1 - e s ) + (t - 0.5) 

Obviously, y = and y = 1 represent pure fluid and pure solid states, respectively. 

The solid volume fraction e s is related to the porosity of lattice cell S2d, namely S2d = 1 _ £s 
. E2d, an important parameter that influences both the gas-phase and the solid-phase motion, 
can be calculated based on the area occupied by the particles in the lattice cell. The porosity of 
two-dimensional plane is given by 

Wtgt ■ 7rr 2 

sid = 1 p — (10) 

where r is the radius of solid particle and h is lattice spacing between nodes. As the drag-force 
correlation is based on 3D systems, we also use a 3D porosity transformed from 2D porosity for 
this purpose. According to (Ouyang and Li, 1999a), 

e 3 d = l--^L(l-£2d) 3/2 (11) 
V^V3 



In order to capture turbulent structures in gas-solid fluidization systems which are usually 
associated with large Reynolds numbers or become turbulent in nature, meso-scale turbulence 
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models are incorporated into the volume-averaged Navier-Stokes equations in DPS (Zhou et al., 
2004a, b; Gui et al. 2008), here LES incorporated into the LBM is used (Krafczyk et al., 2003; 
Yu et al., 2005). The effect of the small sub-grid eddies on the large-scale flow structures are 
modeled through an additional turbulent viscosity v e . 

In Smagorinsky SGS model, v e depends on the strain rate: 

v e = (C s Ax) 2 ||S||, S = ^jlSijSij (12) 

Here C s is Smagorinsky constant and strain rate tensor 5, ; - = (djUi + diU^jl can be obtained 
directly by computing the momentum fluxes (2<y, which are second-order moments of the non- 
equilibrium distribution function (Yu et al., 2005): 

S u = T^Qij, Qij = J] c ki c kj (f k - f k ") , S = -4-, Q = JlQ^j (13) 

where is the kth component of the lattice velocity c, . Based on an eddy relaxation time 
assumption for sub-grid scale stress, the effect of the flow structure at the unresolved scale is 
modeled through an effective collision relaxation time r e , which is the eddy relaxation time with 
respect to the eddy viscosity v e . The eddy viscosity is then incorporated into the LBE by using 
r t = t + r e instead of t in Eq. (6). Accordingly r t 

3 13 1 

T t = -~ v t + - Af = -r- (v + v e ) + -At (14) 

c l 2 c z 2 

where c is the lattice speed given by h/At . From Eq. (13) and Eq. (14), a quadratic equation is 
obtained, which yields 



r, - t+ r e = l -[r+ ^ + 18C;(A.v) 2 e) (15) 



It should be pointed out that Eq. (15) is suitable for both laminar and turbulent flows. On 
one hand, under high gas velocity, Q is large and hence r t » r , which expresses the effect of the 
motion at the unresolved scale. On the other hand, under low gas velocity, Q vanishes and r t ~ r 
, describing the laminar flow. 

2.2. TDHSfor particle dynamics 

Each individual solid particle, treated as a single and point-volume particle, has four prop- 
erties: mass (m), radius (r), position (P) and velocity (v). The motion of solid particle k is 
computed according to Newtons second law of motion as follows: 

m-£ = mg + (F d ), + (F c ), (16) 
at 

where g is the gravitational acceleration, (F^ is the drag force, (F c ) k is the collision force. 

Time-driven hard-sphere model (Hopkins and Louge, 1991) is used for particle-particle col- 
lisions. In TDHS, particles interact as binary and quasi-instantaneous hard-spheres, and both the 
collisions are detected and the particle state are updated at constant time intervals. Namely, when 
two particles are in contact with each other and if the internal product of Pi - P 2 and Vi - V2 is 
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Figure 2: Schematic illustration of two colliding particles 



negative, which means that the particles are moving closer to each other, they will collide (see 
Fig. 2). For the collision between particle 1 and particle 2, the post-collisions velocities can be 
obtained by the following formulae: 

Vi = vi — (Pi - P 2 ) (17) 

m l +m 2 |Pj_p 2 | 2 

(l+e)/wi(Y 1 -Y2)-(Pi-P 2 ) /p D , nsn 

v 2 =v 2 + — - 7 (P1-P2) (18) 

mi + m 2 |Pj - P 2 | 2 

where "*" means post-collision velocities, e is the restitution coefficient that represents the ratio 
of speeds after (post-) and before (pre-) collision, and the velocities in the right-hand side of the 
equations are pre-collision velocities. The walls involved in the simulations are composed of 
infinitely heavy particles with zero velocity, their collision with the fluid particles also stratify 
Eqs. (17) and (18). In the next time step, the particles move to new positions with their new 
velocities. 

It should be pointed out that particle collision is important in simulating gas-solid fluidization, 
especially for dense gas-fluidized beds. In the above time-driven hard-sphere model, only the 
binary-particle collision mechanism and the normal component of contact force are considered, 
which limits the model to low solid holdup conditions (less than 30~40% by volume). For dense 
gas-fluidized beds, both the multi-particle collision mechanism and the tangential component 
of contact force need to be considered, which can be implemented by discrete element method 
(DEM) (Cundall and Strack, 1979). 
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2.3. Inter-phase momentum transfer 

The drag correlations are of great practical importance in numerical models that predict the 
flow behavior of gas-solid fluidized beds of granular material. The calculation of the inter-phase 
momentum transfer coefficient based on the traditional Ergun and Wen & Yu correlations (Gi- 
daspow, 1994; Gidaspow and Jiradilok., 2009) are valid for homogeneous particle suspensions 
only, and are insufficient to capture the heterogeneous structures of gas-solid fluidization. A 
structure-dependent drag based the energy minimization multi-scale model (EMMS) were used 
here, that is (Yang et al., 2003) 



3 gg(l-gg) 



p g |u - v| Cdo • w (s g ) , e g >0.74 



1 150 ^(fr + 1-75(1 " e e )% lu - v| , e g < 0.74 
where according to Yang et al. (2003), the structure-dependent correction factor 
-0.5760+ V^l* , 0.74<e B < 0.82 



CO 



W- 



4( £g -0.07463V +0.0044' 



-0.0101 + - 2^Q2| 0.82 <s„ < 0.97 (20) 

4(e g -0.7789) +0.0040 h 

-31.8295+32.8295e g , e g >0.97 



The drag coefficient Cdo for isolated rigid spherical particle can be calculated by the Schiller 
and Nauman (1935) equation 



+ 0.15fle°' 687 ) , Re p < 1000 



( 0.44, Re p > 1000 

Here the particle Reynolds number is defined as follows: 

= gg P g lu-v|^ P 

Compared with the calculation of the structure-dependent (i.e. EMMS) drag, the particles 
heterogeneously distributed in gas flows should be taken into account for the fluid solver in a 
similar way, the structure-dependent correction factor is also introduced to the LBE, thus the 
modified LBE can be improved as: 

Mx + c,Af, t + At) = fix, t) - 1(1 - y(£s, T t ))(/Xx, f) - f, e Hx, t)) 
+ co(e g )y(s s ,T)Ql 

The force exerted by the particles on the fluid in a computational cell, F p ^ g , can be calculated 
by summing up the momentum transfer that occurs in all discrete directions as 

Fp-g = ^(*g)ri>» Ci (24) 

i=0 

Then the drag force on solid particle k, (F^) k , is calculated according to the following equa- 
tion (Hoomans et al., 1996; Li and Kuipers, 2003): 



( y P )A 



1 -s k 
8 



Here the local inter-phase momentum transfer coefficient the local gas velocity u* and the 
local voidage s k are all calculated by bilinear interpolation using the values of four surrounding 
grid nodes. The total drag force on the particle in a computational cell, F g ^ p , can be obtained 
through a similar computation, 



According to Newtons third law, F g ^ p =F p ^ g . However, due to model inaccuracy numerical 
errors, this is not automatically satisfied in Eqs. 24 and 26. To correct this inconsistency, the 
drag force on solid particle k, (Fd) k , is then adjusted to 



2.4. Computational algorithm and implementation 

The algorithm to implement the method described above is outlined in Fig. 3, which mainly 
includes the following steps: 

1) Input the initial data such as the size and geometry of the simulation domain (for both the 
solid and the fluid phase), the specified boundary conditions (i.e. particle distribution, initial 
flow field, velocity profile, pressure, etc.). The wall boundary condition given in this paper 
is the bounce-back scheme for fluid and specular reflections for solid particles (Wang et al., 
2006). 

2) Perform the statistical calculations of the average velocity and the porosity of solid phase at 
the fluid lattice nodes according to Eq. (8) and Eq. (11), respectively. 

3) Perform the consecutive propagation and collision process over a discrete fluid lattice cell, 
and then evolve the fluid flow according to the modified LBE (Eq. (23)), finally the flow field 
can be obtained from Eq. (5). 

4) The local gas velocities, the local porosity and the local inter-phase momentum transfer co- 
efficient at the position of the center of the particle are calculated by bilinear interpolation 
using the values of the surrounding computational lattices. Solve for F p ^ g from Eq. (24) and 
Fg^ p from Eqs. (25) and (26) at all lattice nodes. Substitute these values of F p ^ g and F g _, p 
into Eq. (27), solve for the new values of fluid drag force acting on the each particle (Fa)*. 

5) Integrate the particle motion equation (Eq. (16)) numerically to calculate the velocities and 
the trajectories of individual particles. 

6) Detect the instantaneous collisions of solid particles according to new particle velocities and 
positions. If the instantaneous collisions occur, the post-collision velocities of individual 
particles are determined according to Eqs. (17) and (18). 

7) Evolve the next time and go back to step 2 unless the set time step is reached. 

3. Simulations and results 

To validate the proposed model, two typical fluidization processes, namely a single bubble 
injected into a fluidized bed at incipient fluidization condition and a particle clustering in the 
riser of a circulating fluidization bed, are simulated and compared to experimental data. 




(26) 




(27) 
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Input initial 
data 



Initial condition (density, velocity and position of both solid 
particles and fluid) and specified boundary conditions 



Calculating porosity and averaged velocity of solid 
particles in a cell 



Collision of distribution function at the lattice nodes 



Propagation of distribution function at the lattice nodes 



f 

Calculating interphase momentum transfer coefficient at 
the lattice nodes 

I 

Calculating local averaged values of porosity, fluid 
velocity, drag coefficient by bilinear interpolation 




Figure 3: The flow chart of LBM-based discrete particle simulation. 
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Table 1: Simulation and physical parameters for 2D bubbling simulation 



W(m) 

H(m) 

nnj(m) 

v inj (m/s) 

v mf (m/s) 

d p (m) 

p g (kg/m 3 ) 

p s (kg/m 3 ) 

/^ g (kg/(ms)) 

f(s) 



Bed width 



0.00675 



e 



Bed height 

Orifice width 

Injection velocity 

Minimum fluidization velocity 

Solid particle diameter 

Gas density 

Solid density 

Gas viscosity 

Time step 

Coefficient of restitution 



0.027 



0.00081 



3.6 



0.189 
5.4 x 10~ 5 



930.0 

1.8872 x 10" 5 
6.0 x 10~ 6 
0.9 

2.7 x 10~ 4 
15000 



1.1795 



h(m) 
N 



Particle number 



Cell or lattice size 



3.1. A single bubble injected into a fluidized bed at incipient fluidization 

In this example, we examine the bubble formation when a pulse of gas is injected at the bot- 
tom of a fluidized bed at incipient fluidization. The data used for the simulation are summarized 
in Table 1 . 

Figure 4 presented the snapshots of the bed at different moments after bubble injection. We 
observed the growth of the spherical cap bubble from an initially small perturbation in the bulk 
of the bed (see Fig. 4(a)~(e)) to its subsequent detachment and rising (see Fig. 4(f)~(h)). When 
the kidney-bubble detached from the bottom of the bed, the bubble wake was followed (see Fig. 
4(j)~(n)). When bubbles arrived at the surface, it lifted some particles from wake into the bubble 
and caused the break-up of the bubble (see Fig. 4(o)~(p)). This phenomenon is in qualitative 
agreement with previous simulation and experimental results of bubbling (Rowe and Yacono, 
1976; Nieuwland et al., 1996). It should be noted that the rising bubble shape is more spherical- 
cap than the elliptic bubble shape in the experimental and numerical results of (Bokkers et al., 
2004) due to the smaller solid particle diameter. 

Figure 5 shows the gas phase leaves the roof of bubble, and the downward moving particles 
near wall drag the gas to the bottom of the bubble which it re-enters the bubble region, which 
results in a pair of symmetrical vortices can be observed in the neighborhood of the rising bubble. 
In additional, the simulated fluidizing gas streamlines round a rising bubble are close to that 
predicted by the bubble model of Davidson and Harrison (1963). These results illustrate the 
modified LBE is the ability of capturing the detailed flow structure of gas phase. The computed 
volume-averaged equivalent bubble diameter De at different injection velocities is shown in Fig. 
6 together with the correlations of Cai et al. (1994). The results show that De increases with the 
increasing of injection velocities, and the deviation is below 5% suggesting that the simulation 
results are reasonable. 
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0.38m/s 



Gas phase velocity: | 
Figure 5: Snapshots of the flow field of gas phase at simulation time t=0.035s. 



q° 0.004 



Simulation (this work) 
- Cai et al.( 1 994) correlations 



_l i l_ 



2.0 3.0 4.0 5.0 6.0 

v.Xm/s) 



7.0 8.0 



Figure 6: Comparison of the correlation formula (Cai et al., 1994) obtained and simulated bubble diameter for a range of 
injection velocities. 
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Table 2: Simulation and physical parameters for 2D riser of circulating fluidized beds 



W(m) 


Bed width 


0.0162 


H(m) 


Bed height 


0.194 


o s 


Initial solids volume fraction 


0.09007 


d p (m) 


Solid particle diameter 


5.4 x 10~ 5 


p g (kg/m 3 ) 


Gas density 


1.1795 


Ps(kg/m 3 ) 


Solid density 


930.0 


yu g (kg/(ms)) 


Gas viscosity 


1.8872 x 10 


f(s) 


Time step 


6.0 x 10~ 6 


e 


Coefficient of restitution 


0.9 


Vinlet(m/S) 


Gas velocity at the inlet 


1.52 


h(m) 


Cell or lattice size 


2.7 x 10~ 4 


N 


Particle number 


123921 



3.2. Riser behavior of fluidization systems 

The riser hydrodynamics of circulating fluidized beds have been investigated extensively 
in both experiments and simulations in past decades. The conditions for this simulation are 
summarized in Table 2. 

Figure 7 shows some snapshots from this simulation. Apparently, radial and axial heteroge- 
neous structures are formed spontaneously and gradually from the homogeneous initial state and 
are present all the times hereafter. 

Figure 8 shows the evolution of the simulated solids flux with time and its comparison with 
experimental data of Li and Kwauk (1994). Obviously, the simulated solids fluxes (13.8kg/(m 2 s)) 
are in a good agreement with the experimental value (14.3kg/(m 2 s)). 

Figure 9 shows the comparison of the axial voidage profiles between the simulated results 
and the experimental data of Li and Kwauk (1994). The simulated profiles are calculated based 
on time-averaged from 4.0s~6.6s. The predicted sigmoid distribution of voidage in the axial 
direction is in reasonable agreement with experimental results, and the deviation may partly be 
ascribed to the unrealistic setup of inlet and outlet boundary conditions in the 2D simulation. 

Figure 10 shows the comparison of the radial voidage profiles between the simulated results 
and the experimental correlations proposed by Tung et al. (1988). The simulated profiles are 
also calculated based on time-averaged from 4.0s~6.6s. The coexistence of a dense annulus and 
a dilute core can be recognized, and the radial profile in the top region is larger than that in the 
bottom region. This means that we get a region with a coexistence of a dilute section at the top 
and a dense section at the bottom. The computed results are in good agreement with empirical 
correlations in the core region, but overpredict the voidage in the annulus region. 

4. Discussion and conclusions 

In this paper, a new simulation model is presented for discrete particle simulation of gas- 
solid fluidization in which the lattice Boltzmann method is applied to model the gas flow, the 
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Figure 7: Evolution of flow structures in the riser of circulating fluidized beds. 
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Figure 8: Evolution of the outlet solid flux in the riser of circulating fluidized beds. 
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Figure 9: Axial voidage profile in the riser of circulating fluidized beds. 




Figure 10: Radial voidage profiles in the riser of circulating fluidized beds. 
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time-driven hard-sphere model is employed to describe the motion of individual particles and 
gas-solid coupling is realized by the EMMS drag. The bubble formation of a single jet into the 
fluidized bed and the characteristic heterogeneous flow structures in the riser of a circulating 
fiuidized bed are simulated successfully by the proposed approach. In general, the agreement 
between the simulated results and the correlations seems to be reasonable. 

From the computational point of view, the physical values of one lattice (i.e. the characteristic 
lattice unit) and the relaxation time should be located in a certain range due to the numerical 
stability of LBM. If air is considered in gas-solid fluidization, the physical length corresponding 
to one lattice may be chosen below the millimeter scale. To guarantee the accuracy of statistical 
calculations for both the average velocity and the porosity of solid phase at the fluid lattice 
nodes, the fluid grid resolution in DPS demands not less than five particle diameters. Therefore, 
the new method is more suitable to simulate the size of particles smaller than millimeter scale, 
such as fluid catalytic cracking (FCC) catalysts for petroleum refiners. On the other hand, from 
the physical point of view, according to hydrodynamics similarity criterion, the new method is 
also capable of simulating the corresponding large actual system only with the requirement of 
matching dimensionless numbers well, such as particle Reynold number, Mach number, Froude 
number, etc. 

In addition, as a linear combination of a direct algorithm and a relaxation scheme, LBM is 
an explicit finite difference approximation to the velocity-discrete Boltzmann equation with a 
collision operator of relaxation type (Junk, 2001). Though LBM for solving fluid flow in the new 
method, in general, is second-order accurate in space and time, it preserves the advantages of ex- 
plicit methods in terms of fast speed and parallelism. On the other hand, time-driven hard-sphere 
model can be viewed as a combination of molecular dynamics (MD) and direct simulation Monte 
Carlo (DSMC) (Bird, 1994), which inherits the virtue of computational efficiency in DSMC and 
physical picture in MD. Although friction and rotation are not considered in TDHS now, it is a 
reasonable treatment for particles in suspension where particle-fluid interactions are dominant. 
Therefore, it can be expected that the proposed approach, combining LBM and TDHS, is a fast 
numerical method, especially for highly parallel computing. 

For simplicity of implementation, the current simulation lies in the restriction of two-dimensional 
fluid and three-dimensional individual spherical particles, and then it can be regarded as a quasi- 
three-dimensional simulation. Further, its extension to three-dimensions is in principle straight- 
forward, namely using the three-dimensional lattice Boltzmann models (i.e. D3Q15, D3Q19 and 
D3Q27) (d'Humieres et al., 2002) for gas flows. For three-dimensional physical problems of 
interest, despite of the relatively high numerical efficiency of LBM and TDHS, it still requires 
excessive computational cost, for which we will naturally resort to parallel computation due to 
the advantage of the natural parallelisms inherent in LBM and TDHS. Therefore, implement- 
ing three-dimensional massive parallel computation of the proposed model is highly desirable in 
future. 

Finally, we would like point out that the proposed LBM based fluid flow solver can be com- 
bined with other solid phase solvers, such as DEM, DSMC, smoothed particle hydrodynamics 
(SPH) (Gingold and Monaghan, 1977; Monaghan, 1992; Xiong et al., 2011) and particle-in-cell 
(PIC) (Snider, 2001) methods. We would like to extend the present method to complex cases 
involving non-spherical particles in the future work. 

Nomenclature 



j3 drag coefficient, kg/(m 3 s) 
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Af 


time step,s 


Ax 


grid size, m 


7 


weighting function 


F 


force, m/s 2 


g 


gravitational acceleration,m/s 2 


P 


position of solid particle, m 


u 


gas velocity, m/s 


V 


particle velocity, m/s 


X 


position of lattice, m 


jUg 


gas viscosity, kg/(ms) 


V 


kinematic viscosity, m 2 /s 


O 


collision term 




lattice weight factor or structure-dependent correction factor 


Pg 


gas density, kg/m 3 


Ps 


solid density, kg/m 3 


C 


unit velocity of lattice, m/s 




initial solids volume fraction 


£ 


voidage 


£2d 


two dimensional porosity 


S3d 


three dimensional porosity 


£ S 


gas volume fractions 


°s 


solid volume fractions 


C 


lattice speed, m/s 


Cdo 


drag coefficient, kg/(m 3 s) 


c s 


Smagorinsky constant number 


c s 


speed of sound,m/s 


d p 


solid particle diameter, m 


e 


coefficient of restitution 


f 


distribution functions 
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particle equilibrium distribution function 


H 


bed height, m 


h 


cell or lattice size, m 


m 


mass of solid particle, kg 


N 


total number of solid particles in the simulation 


n tot 


number of particles in a cell 


On 


momentum fluxes 


r 


radius of solid particle, m 


Hnj 


orifice width, m 


Re v 


particle Reynolds number 




strain rate tensor, s _1 


t 


time,s 


Vinj 


injection velocity, m/s 


Vinlet 


gas velocity at the inlet, m/s 


Vmf 


minimum fluidization velocity, m/s 


v P 


volume of particle,m 3 


w 


bed width, m 


Acronyms 


BGK 


Bhatnagar-Gross-Krook 


CFD 


Computational Fluid Dynamics 


DEM 


Discrete Element Method 


DPS 


Discrete Particle Simulation 


DNS 


Direct Numerical Simulation 


DSMC 


Direct Simulation Monte Carlo 


D2Q9 


Two-Dimensional 9-Velocity Lattice Boltzmann Model 


D3Q15 


Three-Dimensional 15-Velocity Lattice Boltzmann Model 


D3Q19 


Three-Dimensional 19-Velocity Lattice Boltzmann Model 


D3Q27 


Three-Dimensional 27-Velocity Lattice Boltzmann Model 


EMMS 


Energy Minimization Multi-Scale 
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rUU 


Fluid Catalytic Cracking 


i a 
LGA 


Lattice Gas Automata 


1 DC 

Lbt 


Lattice Boltzmann Equation 


1 RM 


Lattice Boltzmann Method 


LES 


Large Eddy Simulation 


MD 


Molecular Dynamics 


PIC 


Particle-In-Cell 


QBMM 


Quadrature-Based Moment Methods 


SGS 


Sub-Grid Scale 


SPH 


Smoothed Particle Hydrodynamics 


TDHS 


Time-Driven Hard-Sphere 


TFM 


Two-Fluid Model 
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